J ; 


T . 


T w's - 


433 


Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 
flCIDES-Un Editor: G.S. Dulikravich. Washington D.C.. October 23-25. I9&L — 


N9 2 - 13 9$0 


AIRFOIL OPTIMIZATION 
WITH EFFICIENT GRADIENT CALCULATIONS 


Thomas Sorensen 

Department of Aeronautics and Astronautics 
Massachusetts Institute of Technology 
77 Massachusetts Avenue 
Room 37-350 
Cambridge, MA 02139 

1 ABSTRACT 




The viscous airfoil design/analysis code XFOIL was extended to allow optimization using con- 
formal mapping coefficients as design variables. The optimization technique employed was the 
Steepest Descent method applied to a Penalty Function. The gradients of the aerodynamic vari- 
ables with respect to the design variables were cheaply calculated as by-products of XFOII/s 
integral boundary layer Newton solver. The speed of the optimization process was further in- 
creased by updating the Newton system boundary layer variables after each optimization step 
using the available gradient information. Two examples are presented. 


2 INTRODUCTION 

Airfoil design can be broken into two schools of thought. The more recent of the two involves 
the use of inverse design methods whereby the airfoil geometry is generated to match a specified 
pressure distribution. The drawback is in determining what makes a good pressure distribution. 
Many examples of inverse design techniques exist in the literature [1, 2, 3]. The older design 
practice uses trial and error geometry guessing. Each new geometry is evaluated using an airfoil 
analysis method and is compared to previous designs. This is continued until an acceptable 
design is iteratively converged upon. This is a time consuming process, but, it does lend itself to 
numerical optimization techniques. Many methods have been tried for inviscid airfoils, several 
examples of which are given by Vanderplaats [4, 5]. Optimization can be computationally 
intensive, so to be a viable design tool the optimization method employed must be efficient. 
Optimization efficiency can be increased by the use of gradient information but calculation 
of this information adds to the computational burden. One method of obtaining the gradient 
information is to perform finite difference calculations, however, this can be extremely expensive. 

The object of the present research was to modify an existing 2D airfoil design/analysis code 
to calculate gradient information during the analysis procedure, with a minimum of excess 
work, such that this information can be used in an optimization process. The optimizer written 
for the design code was simple and robust, but not necessarily the most efficient since the 
emphasis was on developing the ingredients for the optimization: design variables and gradient 
information. The code used was Drela’s XFOIL code [6]. XFOIL has several design routines, 
and includes both viscous and inviscid analysis routines. Principles from both the design and 
viscous analysis routines were combined to allow viscous optimizations. 

The outline for the remainder of this paper is to first present the governing equations, 
the choice of design variables, and how these variables allow efficient gradient calculations. 
These same gradients can also be used to further speed the optimization process which will be 
presented next. Two design examples will be given at the end. 
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^ 3 ANALYSIS 

3.1 Governing Equations 

The optimization scheme utilized in XFOIL was an iterative ’Steepest Descent’-type. In order to 
use this technique the Objective Function and constraints were combined into a Penalty Func- 
tion such that the constrained airfoil optimization problem is converted into an unconstrained 
problem. A constrained airfoil optimization problem can be stated in Penalty Function form as 

Minimize : P(x) = F(x) + i Kj (ffj(x)) 2 , (1) 

1 j=i 

where. 


9j{x) > 0 for j = l,m 
are the constraints that the airfoil is subject to, and 


( 2 ) 


K : 


0 ffj(x) > 0 

* ?;(x) < 0 


(3) 


are the switches that turn the constraints on and off. The cost parameter, k, is a large positive 
quantity used to control the influence of the constraint on the optimization process [5], The 
Objective Function, F(x), is the function that the optimizer will drive to the lowest possible 
value, subject to the stated constraints, using the design variables x. For airfoil optimization 
the Objective Function could be simply the drag coefficient or a combination of several airfoil 
characteristics such as the negative of the range parameter, — MCj/Cj. 


3.2 Design Variables 

The unit circle in the (.'-plane can be mapped to an airfoil in the z-plane bv the transformation 

[3] 


Q r 1 f OO 'j 

^ = (1- T) (1_, " ) exp|^(^ n + iS n )<;- n |, 71 = 0,1,2,--- (4) 

where, xe te is the trailing edge angle. The design variables employed in XFOIL’s optimizer are 
a finite number of the real and imaginary parts of the complex coefficients of Eq. 4: 


x = {A 2 , A 3 , -. A Na , B 2 , B 3 , - - B Nb } t . (5) 

Using the above notation, there are a total of (N A - 2) + ( JV B - 2) design variables. Each design 
variable corresponds to a single design mode such that the optimal airfoil is constructed by a 
sum of these design modes. A particular convenience of these design variables is that the A n ’s 
control the thickness distribution of the airfoil and the B n ' s the camber distribution. Due to 
this distinction the A n ’s and B n 's will be referred to, respectively, as the symmetric modes and 
the anti-symmetric modes. The first 3 symmetric and anti-symmetric design modes are shown 
in Fig. 1. The solid lines for the symmetric modes indicate the airfoil surface for one value of 
A n . The dashed lines show how the surface (i.e. the thickness) changes as another value of 
A n is used. For the anti-symmetric modes, the lines are not the airfoil surface, but the camber 
lines. The first usable design modes are A 2 and B 2 since A 0 , A X ,B 0 , and Bi are constrained by 
Lighthill’s constraints [2] and therefore are not available as design variables. 
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The .4 n and B n coefficients completely control the airfoil geometry with the exception of 
the trailing edge angle and gap. For a typical airfoil only the first twenty or so C n ’s are required 
to define the airfoil. The value of the design variables for a DAE11 airfoil are plotted in Fig. 2 
as an indication of their magnitudes for a typical airfoil. The higher frequency modes quickly 
become unimportant. In both cases, only approximately the first 15 modes are important. 
The DAE11 geometry is shown in Fig. 3 for reference. The higher modes, however, become 
important for airfoils with small leading edge radii. 

3.3 Aerodynamic Quantities 

For optimization efficiency it is imperative that gradient information be calculated and cal- 
culated cheaply. The gradient information will also prove useful in making XFOIL’s viscous 
analysis procedure run faster as will be shown shortly. 

In its unmodified configuration XFOIL solves a viscous flow around an airfoil by constructing 
3 linearized boundary layer (BL) equations at each airfoil and wake node ( N airfoil nodes, N w 
wake nodes) and solving the resulting system using a Newton solver. For a viscous airfoil 
analysis all aerodynamic quantities of interest are functions of the five BL variables: C T » 9, 
m = u r S *, u „, and S'. In this text C T will represent two quantities: in laminar regions it will be 
the amplitude of the most-amplified Tollmien-Schlichting wave, and in turbulent regions it will 
be the maximum shear coefficient. The Newton system only solves for three of these variables, 
C T , 9, and m, since u e and 6* are related to the first three variables. For more details of XFOIL, 
see Drela [6]. 

To calculate the required BL variable gradients, consider the Newton System used in XFOIL 

MW *-{*}■ (6) 

This equation is a block matrix equation where the i t,l -row, J ^-column block of the Jacobian 
Matrix is 



$*- §£■ 

L (7U Tj Vv j OTTXj 

The corresponding z^-row block of the vectors are 


w = < 

] 

• , {*,} = < 

(M 

9i f 

1 

[ <5m, J 

! 

1 ^ J 


Many of the terms in the Jacobian Matrix are zero, but the detailed structure is not important 
here. 

Equation (6) is constructed using 3 BL equations at each node all with the functional form 
Ri — — li ^1* ^2^ * * ’ i + (9) 

where, R x can be g x , or and the subscripts indicate which node is being considered. The 

edge velocity, u e , is composed of an inviscid and a viscous source contribution, 

L 



436 

Third International Conference on Inverse Design Concepts and Optimization in Engineering Sciences 
(ICIDES-Iin. Editor C,.S. Dulilcravich. Washington D.C.. October 23-25. 1991. 

r i 

= Qk + YL dk i m i' ( 10 ) 

j 

where, the tnvisctd part q k depends on the airfoil geometry and hence A n and B n . The mass 
defect, m, therefore also depends on A n and B n , and so does the viscous residual R, in Eq. (9). 
Consequently, a new Newton system is obtained in the form 


t ' 1 ^{ 1 } = -w- 

The i^-row block of the Jacobian addition, [A], is 



' *L. 

dAi d Ay 


£& a 7 • • • 

' 

BB Kg 



l A i] = 

8 Ay dAy 

■ . . dg, 

dA ”A 

day Bay 

zip- 


(12) 


dK dh. 

oT 3 j BAy 

Bh t 

dh, dh. 

bb\ 

Bh x 

^7 J 



The added vector term contains the changes in the 

design variables 




{A} = { A A 

2 , A, 4 3 , 

1 

Ai?2, AB3, • • 

• } r , 

( 13 ) 


where A( ) imphes a change in the design variables between the current optimization step 
and the next optimization step. The modified Jacobian matrix, [7 | A], is no longer square, 
but during normal viscous calculations the geometry is fixed and thus the A A n and Afl n ’s are 
known (i.e they are zero). Therefore, rewriting Eq. (11) with all knowns on the right* hand 
side and then pre-multiplying both sides by [7] 1 the system reduces to 

W = -W 1 W + [i?]{A}, (i 4) 

where, 


[Z7] = - [7]- 1 [A]. (15) 

The viscous solution is obtained when the residual, {R}, is zero. Thus, at convergence 
Eq. ( 14) will have the same form as a first order Taylor series expansion of the 3 BL equations 
in terms of the design variables. For example, the Taylor expansion for C r , 0 , and m at the »** 
node is 


' 6C t% 

*A 

dC Ti 

&An 

n b 

dCr , ) 

60, 

y = E AA * < 

n= 2 

&A n 

' + £ AB n < 

n= 2 


, , 


dm , 




The Taylor coefficients are the BL variable derivatives being sought and after close examination 
it can be seen that they are the columns of [Z?]. For example, the t tfc -row block of \D) is 


L 


j 
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r 

TXT 

9Cr t 

TXT 

dCr x 

5A Na 

dC Tt 

TBi 

dCr t 

dCr t 

A-] = 


gg. 

TI7 * ‘ 

30, 

■ 5^7 



d*i 

™tr B 


dm t 

Jit 

dm t 

TXT ’ 

• wr 1 - 

dA ”A 

li 1 

O Dj 

fgi •• 

OD$ 

■' vtr*- 
dB 


The elements of this matrix are found not by carrying out the matrix multiplication as 
indicated in Eq. (15) but by solving the original Newton system with the columns of [A] added 
as extra right hand sides. Since a direct matrix solver is used, very little extra work is needed 
to calculate the required sensitivities. In addition, the extra right hand sides only have to be 
included after convergence of the system, not every time the system is solved. 

The above derivation presents a scheme to compute the BL variable gradients if the gradients 
of the BL equations, Eqs. (9), are known (i.e. if the terms of [A] are known). The terms in [A] 
are found by use of the chain rule and are included here without derivation 


dR t _ ( dRi \ td qi - 1 \ (dRA / dq % \ 
dA n \dqi-J \ dA n ) + \dq x ) \dA n )' 

where, 


(18) 


dRj _ dRj 

dqi-x 


dR t 

M'-x 


m i - 1 


ut 


(19) 


is found using Eq. ( 10 ) and the definition of the mass defect, m = u e 6 *. Similarly for the B n 
derivatives. In the above four equations can be /», < 7 ,, or h x . At node t the derivatives depend 
only on the information at that node and the upstream node t - 1 . All the terms in Eq. (19) 
are already available once XFOIL constructs the Newton system. Further details of the above 
equations can be found in the author’s Master’s Thesis [ 7 ], 

The only remaining unknown sensitivities in Eq. (18) are the derivatives of q . These can be 
calculated analytically from the expression for q obtained after the complex potential is mapped 
from the circle-plane to the airfoil- plane. At any point, in the circle-plane, the physical speed 
is 


q = exp < 


In 1 


- £)'“ {e~ ia + * iQ C 1 )) - f >» + iB n )C n J 


( 20 ) 


The derivatives of this equation are remarkably easy and cheap to compute: 



3.4 Geometry Gradient 

Now, all aerodynamic variables that depend on the flow solution have been differentiated, and 
only one further piece of gradient information is necessary; the geometry sensitivity. This 
can be found analytically using the integrated form of Eq. (4), however, in practice there is a 
complication. The difficulty arises due to the need for the geometry gradient for the unit chord 
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' 8irf ° n - E 1“* tion ( 4 ). wh«„ integrated, does no, produce a uni. chord airfoil . nd , h „ tfo „ „] 

gradient will not be for a uiut chord. The geometry is subsequently normalized, however this is 
completely satisfactory for the gradient due to movement of the leading edge. This is not a 
concern for symmetric airfoils and is a relatively small effect for cambered airfoils. Therefore, 
he movement of the leading edge point was ignored in calculations for the gradient of z. 

3.5 Updating BL Variables 

Th ' , N 7'°” sy,,,n ! of XF0IL variable, of ,h, previous solution a, the starting 

pom, of the new solution, therefore, the speed of the optimisation tan be increased by simply 
approx, matmg the BL ..ri.ble, of the new airfoil. This can be done b, adding the following 

per urbations to the BL variables at the old optimiration step at those nodes not affected bv 
the transition point: 


(23) 


W = [0]{A}. 

The AA n ’s and AB n ’s in the {A} vector of Eq. (23) are the changes in the design variables 
between the current and new optimization steps, and are calculated from Steepest Descent 
Equation. The remaining two perturbations, 6u e and 66*, can be found using 


N 


and 


n=2 


Na 


N 


<“• = T. + £ 


n= 2 


° du f 
dB n 


66* = y i£ 

t 2 dA n —^ 2 dB 


. . ^ 96* 


(24] 


(25) 


For a reasonable optimization step size this linear extrapolation will give a good approximation 
new BL variables. Thus, the Newton system constructed during the analysis of the new 

coXon COnVCr8e than ^ n ° UP<lating WCTe d ° ne SinCC il wiU have a better ixritial 

Movement of the upper and lower surface transition points from one panel to another will 

tTTrZv ' ""V C ^T ^ thC VariablC8 that this Unear extrapolation will not work near 
transition points. If not considered separately, the poor transition point approximations 

TfTh t 6 en ? Ugh t0 nCgate thC gainS ^ effidency P romised by the updating. The new location 
of the transition points is approximated and then the BL variables at each panel the transition 

pom s have passed over are ’fudged’ . This ’fudging’ process will only affect the rate at which 
the Newton system converges, it will not affect the converged solution. For CV, 9 , and u the 

approximation across the transition point shift is a Unear extrapolation from the previous two 
approximated points, i.e. 


L 


7 * (26) 
where t is a BL node the transition point has passed over. The equations for 9 and u, are similar 

to or ,e! r=r g ind !* L -T bles ;.? r d u was found to *>* a bett «- » 

. * . *. 1 * All that remains to be able to use these transition point 

approximations is to determine how far the transition point has shifted. This is done using 

dZtran C/ ~, dXtran dXf rnn dXf ran 


^Iran — 


dC T 


-6C r + 


d9 


■69 + 


dx_ 


66 * + 


du t 


-6u t . 


(27) 
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All the derivative terms in the above are already calculated in XFOIL to construct the Newton 
system, so the derivation is complete. 

The convergence histories for a simple test case with and without updating the BL variables 
are shown in Fig. 4. The number of iterations for the Newton solver to convergence is plotted 
versus the optimization step number. The amount of time saved is not extensive, but the low 
cost of updating makes it worthwhile. As the optimization continues the savings will be smaller 
since the step sizes are small. 


4 RESULTS 

The two examples presented in this section were run on a DecStation 5000. These examples 
were chosen to show the various properties of XFOIL’s optimizer, they are not designed to be 
realistic design problems. 

5 Example 1 - Cd minimization, M = 0, a = 0° 

The first test case was designed as a simple example to build faith in the optimization code. A 
NACA 0015 airfoil was used as the seed airfoil with Cd used as the Objective Function. The 
only constraint was to keep the angle of attack constant at 0°. The Reynolds Number based 
on the chord was 10 e . The two design variables used were and A 3 . Using only two design 
variables will allow a pictorial representation of the optimization path to be constructed. 

Figure 5 portrays the optimization space for this test case. The contours are of constant 
Cd and a local minimum is located in the upper left corner. The seed airfoil is located out 
of the picture in the lower right corner and the path taken by the optimizer is marked by 
the crosses. Convergence took 24 iterations and approximately 12 minutes. Figure 5 clearly 
shows the larger step sizes in the first five steps, i.e. in the region of large slope. The step 
directions are perpendicular to the contours, as they should be, where the gradients are large. 
As the optimum is neared the step directions start to parallel the contours. This is due to 
the approximations made in the gradient calculations. This is not a detriment since the exact 
mathematical optimum is relatively unimportant. 

From Fig. 6 it is obvious that the largest drag reductions are produced in the first few 
iterations. This is a recurrent observation. Figure 7 compares the optimal airfoil to the seed 
airfoil. Because only two design modes were utilized, the possible change in the airfoil is small. 
However, large changes were made in Cd by modifying the airfoil such that the transition points 
were moved further aft. 

5.1 Example 2 - Cd minimization, M = 0, Cj = 0.5 

The second example optimized the Cd of an airfoil using 7 symmetric and 5 anti-symmetric 
design modes. The seed airfoil was an NACA 3412 and was constrained for a constant lift 
coefficient and a minimum allowed thickness at 95% of the chord. This constraint was necessary 
to prevent negative thickness airfoils. The cost parameter and the Reynolds number were 
k — 100 and Re — 5 x 10 6 . 

This example was stopped after a viscous Newton system was unconverged at the 38 tH 
optimization iteration. The Penalty Function is shown in Fig. 8 . The drag reduction slows 
slightly after 20 iterations but is definitely still headed down when the optimizer was stopped. 
The optimizer was restarted using the last airfoil generated before the Newton system failed as 
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the new seed airfoil. Optimization convergence was achieved after an additional 15 iterations. 
The optimization required approximately 30 minutes. The drag was further lowered from 
C'd = 0.00389 to Cd = 0.00380. The reason for the unconverged Newton system is unexplained 
but it does not invalidate the results of the optimizer. 

The pressure plots of the seed and optimized airfoils are shown in Figs. 9 and 10, respectively. 
The dashed lines in the C p curves are the inviscid solutions and the solid lines the viscous 
solutions. The waviness apparent in the C p curve of the optimized airfoil is due to the fact that 
higher design modes were not used during the optimization. 

Modification of an airfoil design code to use mapping coefficients as the design variables was 
successfully implemented. Gradient information was calculated within the analysis portion of 
the code with a minimum of extra effort. The gradient information was shown to be accurate 
When used in the proper way, the XFOIL optimizer can become a valuable design tool. 
The optimizer should not be used as a ’black box’ to create perfect airfoils but as a designer’s 
tool that will free the designer to become more creative and productive by reducing the time 
spent in iterative design modifications. The ’optimal’ airfoils obtained should be used to give 
the designer ideas for what characteristics the real airfoil should have. 

There were also several areas in which the XFOIL optimizer did not live up to expectations. 
The first is the limited number of design variables that could be utilized. It was found that the 
optimizer should be restricted to N A < 12 and Nb < 12 because the higher mode derivatives 
became inaccurate. This does not allow the generation of completely general airfoils with the 
chosen design variables. This is a disappointment, however the cheap gradient calculations 
made possible by using the mapping coefficients as design variables make up for this deficiency. 
Another disappointment was the temperamental nature of XFOIL’s Viscous Newton solver. 
This does not destroy the promise of the optimizer it only enforces that some care needs to be 
exercised when using the optimizer. 

Another area for future research is the development of design variables that can also control 
the trailing edge angle and gap, and if possible, be completely general. 
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7 NOMENCLATURE 


F 

Objective function 

X 

General design variables 

9j 

Constraints 

m 

Number of constraints 

A n 

XFOIL thickness design variables (symmetric) 

B n 

XFOIL camber design variables (anti-symmetric) 

n a 

Last symmetric design mode used in optimization 

Nb 

Last anti-symmetric design mode used in optimization 

[J] 

Newton system Jacobian matrix 

Ml 

Addition to Jacobian matrix 

{*} 

Newton system unknown vector 

{A} 

Addition to unknown vector 
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Residual vector 

Aerodynamic variables derivative matrix 
Coefficient of lift 
Coefficient of drag 
Mach number 

Reynolds number based on airfoil chord 
Number of airfoil nodes 


/*i 

(TV? rn^ u e , 6 

%tran 

a 

<2 

c = re XUJ 
A 

H ) 

»( ) 

&( ) 


Number of wake nodes 

Node t boundary layer equations 

Boundary layer variables 

Transition point location 

Trailing edge angle parameter 

Angle of attack 

Inviscid surface speed 

Complex circle-plane coordinate 

Difference operator 

Newton system perturbation 

Real part of the quantity in the parenthesis 

Imaginary part of the quantity in the parenthesis 
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Figure 8: Example 2 - Optimization History 
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